nH0 <- 100
nH1 <- 200
phi_parametro <- 0.2
coeficientes <- function(phi) {
# βARMA: the model from Rocha and Cribari-Neto (2009, 2017)
# is obtained by setting coefs$d = 0 and d = FALSE and error.scale = 1 (predictive scale)
list(alpha = 0, beta = NULL, phi = phi, theta = NULL, d = 0, nu = 20)
}# Função para forçar a execução de `BARFIMA.extract` até que ela funcione
so_quero_que_funcione <- function(func, debug = T, max = 1000, ...) {
FINALMENTE <- F
contador <- 0
while (!FINALMENTE) {
contador <- contador + 1
if (debug) print(paste("Tentativa:", contador))
resultado <- tryCatch({
func(...)
}, error = function(err) {
if (contador >= max) {
stop("Deu ruim, número máximo de tentativas atingido")
}
if (any(grepl("\\.btsr\\.extract\\(", deparse(err$trace$call)))) {
# Ignora erro de extração de resíduos
return(NULL)
}
stop(err)
})
if (!is.null(resultado)) {
FINALMENTE <- T
}
}
return(resultado)
}ewma_qcc <- function(amostra_inicial, lambda, amostra_subsequente) {
ew <- ewma(amostra_inicial, lambda = lambda, nsigmas = 3, plot = F, newdata = amostra_subsequente)
registros <- (nH0 + 1):(nH0 + nH1)
list(
ewma = as.numeric(ew$y[registros]),
LI = ew$limits[, 1][registros],
LS = ew$limits[, 2][registros]
)
}gera_amostras_subsequentes <- function() {
lambda <- 0.2
coefs <- coeficientes(phi_parametro)
amostra_controle <- BARFIMA.sim(
n = nH0,
coefs = coefs,
error.scale = 1,
complete = F,
seed = 1337
)
ultima_observacao <- amostra_controle[nH0]
fit <- BARFIMA.fit(yt = amostra_controle, start = list(alpha = 0, nu = 20, phi = 0.1), p = 1, report = F)
phi_estimado <- fit$coefficients["phi"][[1]]
residuos_controle <- BARFIMA.extract(yt = amostra_controle, coefs = coeficientes(phi_estimado))$error
data.frame(
novo_phi = seq(0, 0.8, by = 0.1)
) %>%
rowwise() %>%
mutate(
observacao = list(1:nH1),
dados = list(BARFIMA.sim(
n = nH1,
coefs = coeficientes(novo_phi),
y.start = ultima_observacao,
error.scale = 1,
complete = F,
seed = 1337
)),
residuos = list(
BARFIMA.extract(yt = dados, coefs = coefs)$error
),
qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(ewma < LI | ewma > LS),
total_fora_de_controle = sum(fora_de_controle),
fracao_fora_de_controle = total_fora_de_controle / nH1
)
}
amostras_subsequentes <- so_quero_que_funcione(gera_amostras_subsequentes)## [1] "Tentativa: 1"amostras_subsequentes %>%
group_by(novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
.groups = "drop"
)ggplotly(
amostras_subsequentes %>%
unnest(
cols = c(novo_phi, observacao, ewma, LI, LS)
) %>%
ggplot() +
geom_line(aes(x = observacao, y = LI, color = "Limite Inferior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = LS, color = "Limite Superior"), linetype = "dashed") +
geom_line(aes(x = observacao, y = ewma, color = "EWMA")) +
geom_point(data = . %>% filter(ewma < LI | ewma > LS),
aes(x = observacao, y = ewma, color = "Fora de controle"), size = 2) +
labs(x = "Observação", y = "Caracterísitca", color = "Legenda") +
scale_color_manual(
values = c(
"Limite Inferior" = adjustcolor("red", alpha.f = 0.4),
"Limite Superior" = adjustcolor("red", alpha.f = 0.4),
"EWMA" = "black",
"Fora de controle" = "red"
)
) +
labs(
title = "EWMA para resíduos βAR(1), com λ = 0.2",
) +
theme.base +
theme(
legend.position = "bottom",
legend.title = element_blank()
) +
facet_wrap(~novo_phi)
) %>%
plotly.baseamostras_monte_carlo_gerador <- function(){
numero_de_execucoes <- 100
expand.grid(
k = 1:numero_de_execucoes,
novo_phi = seq(0.2, 0.6, by = 0.1),
lambda = seq(0.1, 0.4, by = 0.1)
) %>%
mutate(
id = row_number()
) %>%
rowwise() %>%
mutate(
amostra_controle = list(BARFIMA.sim(
n = nH0,
coefs = coeficientes(phi_parametro),
error.scale = 1,
complete = F,
seed = id
)),
ultima_observacao = amostra_controle[nH0],
phi_estimado = BARFIMA.fit(
yt = amostra_controle,
start = list(alpha = 0, nu = 20, phi = 0.1),
p = 1,
report = F
)$coefficients["phi"][[1]],
residuos_controle = list(
BARFIMA.extract(
yt = amostra_controle,
coefs = coeficientes(phi_estimado)
)$error
),
dados = list(BARFIMA.sim(
n = nH1,
coefs = coeficientes(novo_phi),
y.start = ultima_observacao,
error.scale = 1,
complete = F,
seed = id + 1337E4
)),
residuos = list(
BARFIMA.extract(
yt = dados,
coefs = coeficientes(phi_estimado),
llk = F,
info = F
)$error
),
qcc = list(ewma_qcc(residuos_controle, lambda, residuos)),
ewma = list(qcc$ewma),
LI = list(qcc$LI),
LS = list(qcc$LS),
fora_de_controle = list(ewma < LI | ewma > LS),
total_fora_de_controle = sum(fora_de_controle, na.rm = TRUE),
fracao_fora_de_controle = total_fora_de_controle / nH1
)
}
if (file.exists("cache_simulacao_ewma.RData")) {
amostras_monte_carlo <- readRDS("cache_simulacao_ewma.RData")
} else {
amostras_monte_carlo <- so_quero_que_funcione(amostras_monte_carlo_gerador)
saveRDS(amostras_monte_carlo, file = "cache_simulacao_ewma.RData")
}amostras_monte_carlo_resumo <- amostras_monte_carlo %>%
group_by(lambda, novo_phi) %>%
summarise(
mean = mean(fracao_fora_de_controle),
min = min(fracao_fora_de_controle),
max = max(fracao_fora_de_controle),
.groups = "drop"
)
amostras_monte_carlo_resumoggplotly(
amostras_monte_carlo_resumo %>%
ggplot(aes(x = novo_phi, y = mean, color = factor(lambda))) +
geom_line() +
geom_point() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
color = "Valores de λ",
title = "Fração de pontos fora de controle por Φ e λ"
) +
theme.base
) %>%
plotly.baseggplotly(
amostras_monte_carlo %>%
ggplot(aes(x = factor(novo_phi), y = fracao_fora_de_controle, fill = factor(lambda))) +
geom_boxplot() +
labs(
x = "Valores de Φ",
y = "Fração de pontos fora de controle",
fill = "Valores de λ",
title = "Fração de pontos fora de controle por Φ e λ"
) +
theme.base +
facet_wrap(~factor(lambda))
) %>%
plotly.base